home *** CD-ROM | disk | FTP | other *** search
/ Micro R&D 5: Mand 2000 / Mand 2000 - Micro R&D CD-ROM Vol 5.iso / bonusstuff / mandtopi.c < prev    next >
C/C++ Source or Header  |  1993-09-12  |  10KB  |  244 lines
  1. /*
  2.     MandToPi.c - Copyright © 1993 CygnusSoftware.
  3.  
  4.     Cygnus Software
  5.     33 University Square, #199
  6.     Madison, WI, 53715
  7.     USA
  8.     CygnusSoft@cup.portal.com
  9.  
  10.     This  program  is  supplied  along  with  Mand2000, the tryware fractal
  11. exploration  program  for the Amiga.  Mand2000 features multi-pass drawing,
  12. animated  zooming,  joystick driving through the Mandelbrot set, full ARexx
  13. support,  1000+  bit  accuracy,  multiple  fractal  equations,  etc.  68000
  14. required.   MandToPi  and  Mand2000  may  not  be  distributed,  but may be
  15. purchased  from  Cygnus Software for $34.95.  Demo versions of Mand2000 are
  16. available on many networks.
  17.  
  18.     MandToPi is based on a discovery made by Dave Boll in 1991.
  19.  
  20.     Boll's rather surprising discovery was that that ubiquitous number, PI,
  21. could  even  be  found in the Mandelbrot set.  Where?  Well, if you look at
  22. the  two  largest bulbs of the Mandelbrot set (the head and body) you might
  23. notice  that  these bulbs are seperated by two spikes or arrows.  One spike
  24. on  either  side.   And you may wonder how close do those two spikes get to
  25. each  other?   How  thin  is  the  `neck'  of  the  Mandelbrot?   Since the
  26. Mandelbrot  set  is  known  to be connected (all those black areas are tied
  27. together  by  microscopic black filaments) the spikes can't actually touch.
  28. But  if  they don't touch, how close do they get?  As it turns out they get
  29. infinitely  close  -  that  is,  there  is just a single mathematical point
  30. separating the two spikes.
  31.  
  32.     The  `real' coordinate of those spikes (also known as the x coordinate)
  33. is -0.75, a nice round number.  What this program does is it tries a series
  34. of  points,  all  with  real  coordinates of -0.75, but each one a tiny bit
  35. closer  to  the  real  axis.   For  each  point it calculates the number of
  36. iterations  it  takes  for  that point to escape.  All of these points will
  37. escape,  because  all  of the points are on the spikes and therefore not in
  38. the Mandelbrot set.
  39.  
  40.     What  Boll's  discovery  was is that there is a very strange pattern to
  41. the  number  of iterations.  Each point we try is one tenth the distance to
  42. the  real  axis,  and each number of iterations that the program prints out
  43. looks more and more like PI!  In fact, the general rule is that if you take
  44. the  point (-.75, E), where E is any small number (less than one), then the
  45. number  of  iterations until that point escapes,, times E, is approximately
  46. equal  to  PI.  The smaller E gets, the closer the approximation.  In fact,
  47. the  error  is  always less than E, so in theory you can get an arbitrarily
  48. accurate  value  for  PI,  although  to get more than seven or eight digits
  49. takes a phenomenal amount of computing time.
  50.  
  51.     Here's  the results from a ten pass run of this program.  This run took
  52. five  hours  on  an  '040  - not surprising as over thirty billion floating
  53. point operations were involved.
  54.  
  55.     E is 1.00000000000, num iterations is           3, PI is ~3.00000000000
  56.     E is 0.10000000000, num iterations is          33, PI is ~3.30000000000
  57.     E is 0.01000000000, num iterations is         315, PI is ~3.15000000000
  58.     E is 0.00100000000, num iterations is        3143, PI is ~3.14300000000
  59.     E is 0.00010000000, num iterations is       31417, PI is ~3.14170000000
  60.     E is 0.00001000000, num iterations is      314160, PI is ~3.14160000000
  61.     E is 0.00000100000, num iterations is     3141593, PI is ~3.14159300000
  62.     E is 0.00000010000, num iterations is    31415927, PI is ~3.14159270000
  63.     E is 0.00000001000, num iterations is   314159266, PI is ~3.14159266000
  64.     E is 0.00000000100, num iterations is  3141592656, PI is ~3.14159265600
  65.  
  66.     The  final  digit  is  sometimes off by a bit, but the other digits are
  67. always  correct.   Strange things these fractals.  Strange stuff math.  You
  68. just can't keep a good number down.
  69.  
  70.     The actual value of PI, to more digits that you'll ever need is:
  71.  
  72.     3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944...
  73.  
  74.     This was NOT calculated with this program.  Don't even try.
  75.  
  76.     For the more graphically oriented, you may want to explore these spikes
  77. with  Mand2000.  However if you aren't careful in your exploration, you may
  78. be  disappointed.   The  first  important  thing to do is to make sure that
  79. Sneaky  pixels  (in  the Setup/Set misc menu) is turned off.  Next, because
  80. the  spikes  are  very  thin  and  can  easily sneak between two columns of
  81. pixels,  you must have the pixels on your monitor lined up so that there is
  82. a  column of pixels at exactly -0.75.  The best way to deal with this is to
  83. zoom  in a bit and then scroll the picture so that the two spikes are quite
  84. close  to  the  left  edge of the screen.  Then bring up the `Set location'
  85. requester  and  type  in -0.75 for the left coordinate.  Now, even when you
  86. scroll  the spikes back to the centre of the screen, or zoom in, there will
  87. still  be  a column of pixels precisely lined up on the spikes.  This makes
  88. it  much  easier to see how close to each other the spikes really get.  The
  89. next step is to increase the number of iterations as high as possible.  The
  90. final  thing  to  do  is to crank up the accuracy a bit - go into the Setup
  91. in  the  Misc menu and set extra accuracy to about four.  Normally Mand2000
  92. tries  to  get away with as little accuracy as possible, for speed reasons,
  93. but  for  calculations  like  this, precise results need a few more digits.
  94. All this should insure that your results will not be severely degraded.
  95.  
  96.     The  source to this program was provided so that mathematical explorers
  97. can  search  around  and  find  other places in the Mandelbrot set where PI
  98. occurs.   It  can  also  be  found,  for  instance,  at  the  `butt' of the
  99. Mandelbrot  set,  but  I'll  leave  the exact method as an exercise for the
  100. reader.
  101.  
  102.  
  103.      To  find  out  more  about  Boll's  discovery,  read  "Fractals for the
  104. Classroom - Part Two", an excellent book for anybody interested in learning
  105. wierd and wonderful facts about these beatiful monsters.
  106.  
  107.     When  compiling  this  program  be sure to compile with a very accurate
  108. floating  point  option.   The `Fast Floating Point' (FFP) routines are not
  109. accurate  enough and will cause the program to run forever.  If you have an
  110. FPU, compile and link with the FPU options, because they are both very fast
  111. and very accurate.  And be sure to link with the appropriate floating point
  112. library, needed to print out the results.  */
  113.  
  114.  
  115.  
  116. #include <stdio.h>
  117.  
  118. #define    AMIGA
  119.  
  120. #ifdef    AMIGA
  121. /* The special code is necessary for the Amiga to allow this program to */
  122. /* be run from the workbench, where a console window must be opened and */
  123. /* Ctrl-C won't stop the program. */
  124. #include <clib/exec_protos.h>
  125. #include <dos/dosextens.h>
  126. #endif
  127.  
  128. /* This may be an Aztec specific function - replace it with your own or */
  129. /* just comment it out.  It just checks for control C and exits when it */
  130. /* detects it. */
  131. long Chk_Abort(void);
  132.  
  133.  
  134.  
  135.  
  136. /*
  137.     Warning.   Set  this  to higher numbers at your own risk.  Even at just
  138. seven,  the  final  pass  will  do 3,141,593 iterations - which takes a few
  139. seconds  even  on  an  '040.   Each  time  you increment NUMPASSES, it does
  140. ten  times  as  many  iterations - if you set it to ten you can expect your
  141. computer to be busy for several hours at least.
  142.     */
  143.  
  144. #define    MAXPASSES    10
  145.  
  146. /*
  147.     For those hard core experimenters who want to try multi-week runs, I've
  148. made  counter  into  a long double instead of a mere long.  This will allow
  149. maximum  counts of somewhere around 16 billion billion instead of only four
  150. billion.   (Half  of  you  will probably complain that a long double can go
  151. much  higher  than  that  -  but  it  can't if you still want to be able to
  152. accurately  add  one  to  it).   However, after that many iterations, it is
  153. unlikely  that  enough  accuracy  will  be left to continue giving you more
  154. digits of PI.  That's chaos for you.
  155.     */
  156.  
  157.  
  158.  
  159. int main(int argc, char *argv[])
  160.     {
  161.     register long double    zr, zi, cr, ci;
  162.     register long double    rsquared, isquared;
  163.     register long double    counter;
  164.     long double starti;
  165.     register short    abortcounter;
  166.     short    pass;
  167.     FILE    *mystdout;
  168.     unsigned short    breakable;
  169. #ifdef    AMIGA
  170.     struct Process *myprocess;
  171.  
  172.     myprocess = (struct Process *)FindTask(NULL);
  173.     if (!(myprocess->pr_ConsoleTask))
  174.         {
  175.         mystdout = fopen("con:0/0/640/150/Pi from the Mandelbrot set./Wait/Auto/Close", "w");
  176.         breakable = FALSE;
  177.         }
  178.     else
  179. #endif
  180.         {
  181.         mystdout = stdout;
  182.         breakable = TRUE;
  183.         }
  184.  
  185.  
  186.     fprintf(mystdout, "MandToPi: Copyright © 1993 Cygnus Software.\n");
  187.     fprintf(mystdout, "Calculating PI based on escape times from the `neck' at (-0.75, 0.0) of\n");
  188.     fprintf(mystdout, "the Mandelbrot set.  Read the file 'MandToPI.c' for more information.\n");
  189.     fprintf(mystdout, "\n");
  190.     if (breakable)
  191.         {
  192.         fprintf(mystdout, "Type Ctrl-C to stop.\n");
  193.         fprintf(mystdout, "\n");
  194.         }
  195.  
  196.     /* The initial value for starti and the amount to reduce it by are quite */
  197.     /* arbitrary.  You could just as easily start at 0.683 and multiply by */
  198.     /* 0.5 each time instead.  PI is unavoidable. */
  199.     for (starti = 1.0, pass = 0;  pass < MAXPASSES && (breakable || pass < 5);  starti *= 0.1, pass++)
  200.         {
  201.         zr = 0.0;
  202.         zi = 0.0;
  203.         cr = -.75;
  204.         ci = starti;
  205.         rsquared = zr * zr;
  206.         isquared = zi * zi;
  207.  
  208.         /* The guts of the program - standard Mandelbrot iteration formula. */
  209.         for (counter = 0; rsquared + isquared <= 4.0; counter++)
  210.             {
  211.             zi = zr * zi * 2;
  212.             zi += ci;
  213.  
  214.             zr = rsquared - isquared;
  215.             zr += cr;
  216.         
  217.             rsquared = zr * zr;
  218.             isquared = zi * zi;
  219.  
  220.             if (breakable)
  221.                 if (!((++abortcounter) & 0xFF))
  222.                     Chk_Abort();
  223.             }
  224.  
  225.         /* To make the results look pretty, be sure to put .N, where is N is MAXPASSES or greater */
  226.         /* after each percent sign in the following printf control string.  Otherwise you may get */
  227.         /* misaligned printouts or, worse, missing digits. */
  228.         fprintf(mystdout, "E is %.11Lf, num iterations is %11.Lf, PI is ~%.11Lf\n", starti, counter, starti * counter);
  229.         }
  230.  
  231.     fprintf(mystdout, "\nDone.\n");
  232.  
  233. #ifdef    AMIGA
  234.     if (!(myprocess->pr_ConsoleTask))
  235.         {
  236.         fprintf(mystdout, "Run from the CLI to see a longer run.\n");
  237.         fclose(mystdout);
  238.         }
  239. #endif
  240.  
  241.     return 0;
  242.     }
  243.  
  244.